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ABSTRACT 

We introduce an exact Bayesian approach to search for non-Gaussianity of 
local type in Cosmic Microwave Background (CMB) radiation data. Using sim- 
ulated CMB temperature maps, the newly developed technique is compared 
against the conventional frequentist bispectrum estimator. Starting from the 
joint probability distribution, we obtain analytic expressions for the conditional 
probabilities of the primordial perturbations given the data, and for the level of 
non-Gaussianity, /nl, given the data and the perturbations. We propose Hamil- 
tonian Monte Carlo sampling as a means to derive realizations of the primordial 
fluctuations from which we in turn sample /nl- Although being computationally 
expensive, this approach allows us to exactly construct the full target posterior 
probability distribution. When compared to the frequentist estimator, applying 
the Bayesian method to Gaussian GMB maps provides consistent results. For 
the analysis of non-Gaussian maps, however, the error bars on /nl do not show 
excess variance within the Bayesian framework. This finding is of particular rel- 
evance in the light of upcoming high precision CMB measurements obtained by 
the Planck satellite mission. 

Subject headings: cosmic background radiation — cosmological parameters — 
methods: data analysis — methods: numerical — methods: statistical 
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Introduction 



Precise measurements of tlie cosmic microwave background (CMB) radiation have vastly 
improved our understanding of cosmology and played a crucial role in constraining the set 



of fundamental cosmological parameters (Spergel et al. 2003, 2007 Hinshaw et al. 2009 



Larson et al. 2010). This success is based on a tight connection between the temperature 



fluctuations we observe today and the physical processes taking place in the early universe. 
Inflation is currently the favored theory predicting the shape of primordial perturbations 



(Guth 



1981 



Albrecht & Steinhardt 



1982 



Linde 



1982 



Starobinskii 



1982). In its simplest 



form, it is driven by a single scalar field in ground state with quadratic kinetic term that 
rolled down a fiat potential slowly. This configuration leads to very small non-Gaussianities 
(see Acquaviva et al. 2003 Maldacena 2003 for a first order, and Pitrou et al. 2010 for 



the full second order calculation). Hence, a clear detection of an excess of primordial non- 
Gaussianity would allow us to rule out the simplest models. Together with constraints on 
the scalar spectral index ng and the search for primordial gravitational waves, the test for 
non-Gaussianity therefore becomes another important means to probe the physical processes 
of the early universe. 

In this paper, we focus on non-Gaussianity of local type, where the amplitude of non- 
Gaussianity is measured by a single parameter, /nl (Salopek & Bond 1990). A common 



strategy for estimating /nl is to evaluate the bispectrum of the CMB (Komatsu et al. 2002 



2003 Spergel et al. 2007; Yadav & Wandelt 2008 Smith et al. 2009). This is usually done 



indirectly via a cubic combination of filtered CMB maps reconstructing the primordial per- 



turbations (Komatsu et al. 2005 Yadav et al. 2007, 2008). This approach takes advantage 



of the specific signatures produced by primordial non-Gaussianity, resulting in a computa- 
tionally efficient algorithm. A variant of this estimator has been successfully applied to the 
7-year data release of the Wilkinson Microwave Anisotropy Probe (WMAP), resulting in 
-10 < /nl < 74 at 95 % confidence level (iKomatsu et al.||2010l). 



The bispectrum estimator used in previous analyses has been shown to be optimal, 
i.e. it satisfies the Cramer- Rao bound (Babich 2005). However, this turns out to be true 



only in the limit of vanishing non-Gaussianity (Creminelli et al. 2007). For a significant 



detection of /nl, the estimator suffers from excess variance, a finding that has also been 
verified numerically ( Liguori et al.||2007D . For the simplified case of a fiat sky approximation. 



neglected transfer functions and instrumental noise, Creminelli et al. (2007) showed that it 



should be possible to construct an improved version of the estimator that is equivalent to a 
full likelihood analysis up to terms of the order (9(l/ln A^pix). 



Bayesian methods for the analysis of various aspects of CMB data have been successfully 



3 



developed in the past, e.g., for an exact power spectrum determination using Gibbs sampling 



(e.g. 


Jewell et al.|2004 


Wandelt et al.|2004 Larson et al.|2007 Jewell et al.|200 


9 ) , to separate 


foreground contributions from the CMB anisotropies (e.g. 


Hobson et al.||1998 


Barreiro et al. 


2004 


Eriksen et al. 


2006, I2008a"b| [Dickinson et al. 200£ 


), or to probe for non-Gaussian 


features (e.g. Rocha et al. 2001j|Efstathiou et al. 2009 EnBlin et al. 2009 Vielva & Sanz 2009). 



They offer a natural way to marginalize over uncertainties e.g. attributed to foreground 
contamination or instrumental effects. This is of particular importance for a reliable analysis 
of weak signals and an advantage over frequentist methods, where no such procedures exist. 



Here, we advance the exact scheme introduced in Eisner et al. (2010) to infer the level of 



non-Gaussianity from realistic CMB data within a Bayesian approach. 

We use simulated Gaussian and non-Gaussian CMB temperature maps to compare 
and contrast the conventional frequentist (bispectrum) estimator with the exact Bayesian 
approach. We show that the latter method does not suffer from excess variance for non-zero 
/nl, and can deal with partial sky coverage and anisotropic noise properties, a feature of 
particular importance for local non-Gaussianity and for any realistic experiment. 

The paper is organized as follows. In Sect.|2| we briefly outline the theoretical model used 
to describe primordial non-Gaussianity. We review the conventional frequentist bispectrum 
estimator and present our exact Bayesian approach to infer the amplitude of non-Gaussianity 
in Sect. |3] Then, we use simulated maps to compare the performance of the newly devel- 
oped technique to the traditional estimator (Sect. |4]). We demonstrate the capability of the 
Bayesian scheme to deal with realistic CMB experiments in Sect. [5j Finally, we summarize 
our results in Sect. |6l 

Throughout the paper, we assume the WMAP5+BA0-I-SNALL cosmological parame- 
ters dKomatsu et al.|2009[ ): Qa = 0.721, = 0.1143, Qb = 0.02256, A2j(0.002 Mpc"^) = 
2.457 ■ 10"^ h = 0.701, = 0.96, and r = 0.084. 



2. Model of non-Gaussianity 

The multipole coefficients aim of the CMB temperature anisotropies are related to the 
primordial fluctuations, 

'^^m = ~ y" ^^^^ ^'^^^ j ^("'' ^)^em 9i{k) ji{kr) + ritm 

= M<l>e^ + nim, (1) 

where is the spherical harmonic transform of the primordial adiabatic perturbations at 
comoving distance r, gi the transfer function in momentum space, and je the spherical Bessel 
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function of order ^. Additive noise is taken into account by n^^, for a compact notation we 
will use the operator M as a shorthand for the radial integral in what follows. Traces of 
non-Gaussianity in the primordial fluctuations will be transferred to the multipole moments 
aim according to Eq. [T} potentially making them accessible to CMB experiments. 

We focus on non-Gaussianity of local type, which is realized to very good approximation 



in multi- field inflationary models as described by the curvaton model (Moroi & Takahashi 



2001 Lyth et al. 2003), or in ekpyrotic/cyclic universe models (Khoury et al. 2001 Enqvist 



& Sloth 2002; Steinhardt & Turok 2002). Here, we can parametrize the non-Gaussianity of 



$ via a quadratic dependency on a Gaussian auxiliary field <I>l, that is local in real space. 



of the form (Salopek & Bond 199^ Gangui et al. 1994) 



$L(r) + /NL[$^( 



(2) 



where /nl is a dimensionless measure of the amplitude of non-Gaussianity and we truncate 
the expansion at third order in $l. 

The Bayesian method presented in the following section takes advantage of the simple 
form of Eq. |2j which links the properties of the primordial perturbations $ to that of a 
Gaussian random field $l- As a result, it cannot easily be generalized to the analysis of other 
types of non-Gaussianity, where no such relation exists. Though this poses an important 
limitation of the method, improved statistical means for the search for non-Gaussianity of 
local type are of particular relevance as the conventional bispectrum estimator is known to 
suffer from large excess variance here. Finally, we explicitly stress the interesting possibility 
to include the cubic term $l the perturbational expansion (Eq. [2]) to obtain simultaneously 
constraints to the next order non-Gaussianity parameter, commonly referred to as (/nl- 



3. Analysis techniques 
3.1. Prequentist estimator 



In the following, we briefly review the fast estimator as proposed by Komatsu et al. 



(2005). This estimator is optimal for uniform observation of the full sky. More general least- 



square cubic estimators have been found for data with partial sky coverage and anisotropic 



noise (Creminelli et al. 2006, see also the review of, e.g., Yadav & Wandelt 2010). 



5, 



prim 



To estimate the non-Gaussianity of a CMB temperature map, one constructs the statistic 
out of a cubic combination of the data, 



■'prim 



drr^ / d'^n A(r,n) 



.r,n] 



(3) 
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The spatial integral runs over two filtered maps, 

A{r,n) = ai{r) ai„, Yi^in) 



B{r,n) = 1 (3e{r) aim >^<?m(n) 



that are constructed using the auxiliary functions 

tt£(r) = - I dkk^ gi{k) je{kr) , 



dkk^V{k) geik)jeikr) 



(4) 
(5) 

(6) 
(7) 



and the inverse of the CMB plus noise power spectrum, C^^. The power spectrum of the 
primordial perturbations is denoted by V{k). Now, we can calculate the estimated value of 
/nl from the statistics Sprim by applying a suitable normalization, 

-1 



fi 



NL 



E 



1 



( ryprim\2 

\^ Jhhh '-ii 



prim 1 



Ul<£2<^3 

where A^^^j^g = 6, when (-1 = (-2 = ^3, 2, when ii = £2 ^3 ^1 7^ ^2 
The theoretical bispectrum for /nl = 1, ^e^i^^y is given by 



and 1 otherwise. 



Kms = 2 ^^.^.^3 j drr' (r)/3,,(r)a,3(r) + /3,3(r)/3,,(r)a,,(r) + /3,,(r)/3,3(r)a,,(r)] , (9) 
where a combinatorial prefactor is defined as 



(2£i + l)(2£2 + l)(24 + l) / ii i2 is 





(10) 



Recently, the Bayesian counterpart of the fast estimator has been developed within the 
framework of information field theory by expanding the logarithm of the posterior probability 



to second order in /nl (EnBlin et al. 2009). Here, the equivalent of the normalization factor 



in Eq. [8] becomes data dependent, accounting for the fact that the ability to constrain /nl 
varies from data set to data set. We will go beyond this level of accuracy and present an 
exact Bayesian scheme in the next section. 



3.2. Exact Bayesian inference 



We now introduce a Bayesian method that, in contrast to the bispectrum estimator, 
includes information from all correlation orders. Our aim is to construct the posterior 



- 6 - 



distribution of the amplitude of non-Gaussianities given the data, -P(/nlM)- To this end, we 
subsume the remaining set of cosmological parameters to a vector 9 and rewrite the joint 
distribution as 



p{d, $L, /nl, 9) = p(rf|$L, /nl, e) p($l|0) p(/nl) P{e) . 



Substituting the noise vector in terms of data and signal, we can use Eq. [T]et seq. to express 
the probability for data d given <I>l, /nl, and 6 up to an overall prefactor 



P(rf|<l>L,/NL,^) oce 



^l/2[d-J\/(<I.L+/NLK-K)))l^^"M'i-M(<I'L+/NLK-(-I'?,»)] 



(12) 



where we introduced the noise covariance matrix A^. The prior probability P($l|^) can be 
expressed as multivariate Gaussian distribution by construction, thus, we eventually obtain 

P{d, <I>L, /nl, 9) cx exp I -1/2 [d - M(<I>L + hd'^l - ^ 

X [rf-M($L + /NL($^-(<f^)))] -1/2 4^^-^*l-/^l/2^L} (13) 

as an exact expression for the joint distribution up to a normalization factor, assuming a 
Gaussian prior for /nl with zero mean and variance cxjj^^ , and a flat prior for the cosmological 
parameters. The covariance matrix P$ is constrained by the primordial power spectrum 



predicted by inflation, V{k), and given by (Liguori et al. 2003) 



$L i.mM) i^m.ir,)) = - 5% 5Z / dk k'Vik) je,{kn) je.ikr. 



(14) 



To evaluate the joint distribution (Eq.[l3| directly would require to perform a numerical 
integration over a high dimensional parameter space. For realistic data sets this turns out 
to be impossible computationally. We pursue a different approach here. First, we note that 



the exponent in Eq. 13 is quadratic in /nl and hence the conditional density P(/nlM, '^'l, ' 



is Gaussian with mean and variance 



(/nl) = ((/nl - (/nl))')($^ - {<^l)yM^N-\d - M<|.l) 
((/nl - (/nl))^) = m - {^l))^M^N-^M{^l - ($i)) + llal,y' 



(15) 



Thus, for any realization of $l, Eqs. [15] permit us to calculate the distribution of /nl given 
the data. Similarly, we can calculate the conditional probability P($lM, ^) by analytically 
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marginalizing Eq. [13] over /] 



INL, 



c(fNLi^($L,/NLM,l 



-l/2(d-M<I>L)t 



X e 



N- 



'lAf(<I.2-(4>2))(*2_ 



(#^))tAftjv- 



x(d-M$L)-l/2 4p-icI.L 



(16) 



Now we can outline our approach to infer the level of non-Gaussianity from CMB data 



iteratively. First, for given data d, we draw "I'l from the distribution Eq. 16 Then, /nl can 



be sampled according to Eqs. [15] using the value of •I'l derived in the preceding step. If the 
samphng scheme is iterated for a sufficient amount of cycles, the derived set of /nl values 
resembles an unbiased representation of the posterior distribution P(/nlM, 0). 

Unfortunately, there exists no known way to draw uncorrelated samples of $l from 
its non-Gaussian distribution function directly. Here, we propose Hamiltonian Monte Carlo 
(HMC) sampling to obtain correlated realizations of the primordial perturbations. Contrary 
to conventional Metropolis-Hastings algorithms, it avoids random walk behavior in order to 
increase the acceptance rate of the newly proposed sample. This is a mandatory require- 
ment to explore successfully high-dimensional parameter spaces as found here. For HMC 
sampling, the variable is regarded as the spatial coordinate of a particle moving in a potential 



well described by the probability distribution function to evaluate (Duane et al. 1987). A 



generalized mass matrix W and momentum variables p are assigned to the system to define 
its Hamiltonian 

(17) 



H = 1/2 p'^W-^p - log[P($LM, 



where the potential is related to the posterior distribution as defined in Eq. 16 The system 



is evolved deterministically from a starting point according to the Hamilton's equations of 
motion 



d$L 
dt 
dp 
dt 



dH 

dp 
dH 



d\og[Pi^L\d,i 
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which are integrated by means of the second order leapfrog scheme with step size 5t, 



p{t + 



6t. 



Pit) + 



St aiog[p($LM, 



9$T 



-I>L(t) 



$L(t + 6t) = <l'L(t) + 6t W-^p{t + 



6t. 



p{t + 6t)=p{t+-) + - 



(19) 



*L(t+5t) 

The equation of motion for $l can easily be solved, as it only depends on the momentum 
variable. To integrate the evolution equation for p, we derive 



glog[P($LM,l 
5$T 



N' 



a 



/nl 



($2))tMtiV-iM($2 - ($2)) + 
X (c/-M$l) 



(20) 



as an approximate expression neglecting higher order terms in $l- The final point of the 
trajectory is accepted with probability a = min(l, exp[— Aif]), where AH is the difference 
in energy between the end- and starting point. As the energy is conserved in a system with 
time-independent Hamiltonian, the acceptance rate in case of an exact integration of the 
equations of motion would be unity, irrespecticive of the complexity of the problem. Intro- 
ducing the accept/reject step restores exactness also in realistic applications as it ehminates 



the error originating from approximating the gradient in Eq. 20 and from the numerical 
integration scheme. In general, only accurate integrations where AH is close to zero re- 
sult in high acceptance rates. This can usually be archived by choosing small time steps 
or an accurate numerical integration scheme. However, as the time integration requires the 
calculation of spherical harmonic transforms with inherently limited precision, higher order 
methods turn out to be unrewarding. Furthermore, the efficiency of a HMC sampler is 



sensitive to the choice of the mass matrix W. In agreement with Taylor et al. (2008), we 



found best performance when choosing W as inverse of the posterior covariance matrix of 
the primordial perturbations, which we derive from the Wiener filter equation for purely 
Gaussian perturbations to good approximation, 

P(rf,<l>^,^) oc exp|-l/2 [c/-M$^]"^A^~^ [rf-M$^] - l/2$^tp^i$^} , (21) 
with mean and variance of the distribution P($'^|(i, 9) 

_ ($G^)2^ _ [M^N-^M - P^^y^ , hence 



p-i 



(22) 
(23) 
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For the calculation of the mass matrix W in the presence of anisotropic noise or partial 
sky coverage, we still adopt a simple power spectrum as approximation for A^~^ in spherical 
harmonic space at the cost of a reduced sampling efficiency. 

We initialize the algorithm by performing one draw of the primordial perturbations from 
the Gaussian posterior P($'^|rf, 9) (Eqs. [22|). 



4. Scheme comparison 

We use simulated CMB temperature maps obtained with the algorithm described in 



Eisner & Wandelt (2009) to compare the newly developed Bayesian scheme to the conven- 
tional frequentist approach. We chose a Gaussian (/nl = 0) and a non-Gaussian (/nl = 100) 
CMB realization at a HEALPix resolution of riside = 256 and imax = 512, superimposed by 
isotropic noise with a constant power spectrum amplitude of = 10~''mK^. We show the 
non-Gaussian temperature map besides the input signal and noise power spectra in Fig. [T} 

Performing the analysis within the frequentist framework, we derive /nl = 4 for the 
Gaussian and /nl = 97 for the non-Gaussian simulation. To obtain an estimate of the 
attributed error, we conducted 1000 Monte Carlo simulations with the input parameters as 
quoted above. For the Gaussian realization, we find a standard deviation of cr^^ = 15, in 
perfect agreement with the value predicted form a fisher information matrix forecast. For the 
non-Gaussian simulation, however, the derived error o"^^ = 20 is already considerably larger 
than in the Gaussian case — the sub-optimality of the bispectrum estimator at non-zero /nl 
becomes manifest. 

In the Bayesian analysis, we construct the full posterior distribution out of the samples 
drawn from it. We chose a Gaussian prior for /nl with zero mean and a very large width 
of o''f"°^ = 500 in order to not introduce any bias to the results. For an efficient sampling 
process, we tuned the time step size 6t of the HMC algorithm to realize a mean acceptance 
rate of about 40 %. To reduce the overall wall clock time needed for the analysis of one CMB 
map, we ran 32 chains in parallel and eventually combine all the samples. For reliable results, 
it is imperative to quantitatively assess the convergence of the Monte Carlo process. Here, 



we apply the statistics of Gelman & Rubin ( 1992 ) to the obtained samples. It compares the 



variance among different chains with the variance within a chain and returns a number in 
the range of < i? < cxd which refiects the quality of the convergence of the chains with a 
given length. In general, a value close to i? = 1 refiects good convergence. As this value 
refers to the convergence of a single chain, we in fact obtain a significantly better result after 
a combination of all of the 32 independent chains we generated. 
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For the Gaussian simulation, we run chains with a length of 25 000 samples each, dis- 
carding the first 5000 samples during burn-in. With these parameters, we find excellent 
convergence as confirmed by the Gelman- Rubin statistics, R = 1.04. The final result along 
with a comparison to the frequentist scheme is shown in Fig. [2j In the Bayesian analysis, we 
find a mean value of (/nl) = 3 and a width of the distribution af^^ = 15. As the bispectrum 
estimator is known to be optimal in the limit of vanishing non-Gaussianity, the two different 
approaches lead to consistent results. 

To repeat the analysis of the non-Gaussian map, we again generated 32 independent 
chains with a length of 40 000 samples each. After dropping the first 10 000 elements to 
account for the period of burn-in, we estimated the convergence of the individual chains by 
means of the Gelman-Rubin statistics and find R = 1.4. The inferred mean of (/nl) = 99 
at an 1-a error of ct/j^j^ = 15 is in good agreement with the input value of the simulation. 
We directly compare the Bayesian to the frequentist result in Fig. [3} where we now find an 
important difference in the outcomes. Whereas for a significant detection of non-Gaussianity 
the frequentist estimator suffers from excess variance, the Bayesian scheme still provides the 
same error bars as for the Gaussian simulation. This increase in variance has been found to 
be an intrinsic property of the conventional bispectrum estimator applied to the detection of 



local non-Gaussianity. Creminelli et al. (2007) show the existence of an improved cubic esti- 



mator which better approximates the maximum likelihood estimator even for non-vanishing 
values of /nl- While this estimator has not yet been constructed for realistic data sets, the 
Bayesian analysis we present here yields as a by-product the maximum a posteriori estima- 
tor which becomes the maximum likelihood estimator in the limit of large prior variance for 
/nl- In addition, the Bayesian analysis produces the full posterior distribution using all the 
information about /nl contained in the data. As we demonstrate in this paper, the variance 
of the posterior distribution does not change in the case of non-zero /nl, but its shape does. 

We note that the computational cost for the Bayesian analysis with the exact marginal- 
ization of the high-dimensional $ parameter space is quite demanding. With the setup as 
described here, the runtime for the Gaussian and the non-Gaussian simulation amounts to 
about 80 000 CPUh and 150 000 CPUh, respectively. It is dominated by spherical harmonic 
transforms that show a scaling behavior of 0{Np(^), where A^pix are the number of pixels in 
the data map. Though computationally expensive, the algorithm in its present implemen- 
tation enables the analysis of WMAP data with an only moderately higher resolution than 
that of the simulations considered here. The reason for the inefficiency of the algorithm lies 
in the large correlation length of the /nl sampling chains. We illustrate this fact in Fig. |4| 
where we display three out of the 32 chains of the non-Gaussian simulation. In addition, we 
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show the autocorrelation function of a chain as defined via 

^(^^) = ^ E ^^^"'^^[.{^"'"^''"^^ . (24) 

i 

where N is the length of the /nl chains with mean /i and variance cx^. 

It is interesting to note that the derived values of /nl and their error bars will in general 
not agree exactly between the two approaches, even for a Gaussian data set. The frequentist 
estimator is unbiased with respect to all possible realizations of signal and noise. The error 
bars, calculated via Monte Carlo simulations, are the same for all data sets with identical 
input parameters by definition. The Bayesian approach, on the other hand, returns the 
entire information contained about the local model in the particular realization subject to 
the analysis. Thus, the uncertainty in the parameter is computed from the data itself and 
will vary from data set to data set, as cosmic variance or accidental alignments between signal 
and noise may impact the ability to constrain the level of non-Gaussianity. Furthermore, 
the Bayesian method constructs the full posterior probability function instead of simply 
providing an estimate of the error under the implicit assumption of a Gaussian distribution. 



5. Application to more realistic simulations 

In the previous section, we have demonstrated the Bayesian approach under idealized 
conditions such as isotropic noise properties and a full sky analysis. However, applying the 
method to a realistic CMB data set requires the ability to deal with spatially varying noise 
properties and partial sky coverage. 

In this context, a general problem is the mixture of preferred basis representations. 
Whereas the covariance matrix of the primordial perturbations can naturally be expressed 
in spherical harmonic space, the noise covariance matrix and the sky mask are defined 
best in pixel space. For the frequentist estimator, this is known to be problematic as e.g. 
in the calculation of the auxiliary map B(r,h) in Eq. [s] (the Wiener filtered primordial 



fluctuations, see also Eqs. 22 for an equivalent, but more didactic expression), the inversion 
of a combination of the two covariance matrices has to be computed. For anisotropic noise, 
this can only be done by means of iterative solvers, whose numerical efficiencies depend 
crucially on the ability to identify powerful preconditioner^ 

For the Bayesian analysis scheme as presented here, however, the relevant equations do 



^This can be very difBcult, see, e.g., the discussion in 



Smith et al. 



(20071 
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not contain any terms of this structure. Therefore, the computations remain straightforward 
even in the presence of arbitrary anisotropic noise properties and sky cuts. To demonstrate 
this abihty, we performed a reanalysis of the simulated non-Gaussian temperature map of 
Sect. |4| now superimposed by anisotropic noise as typically expected for a high frequency 
WMAP channel. With these parameters, the average noise power spectrum roughly remains 
at a level of about C1°'''' ^ IQ-'^mK^ but the noise is no longer spatially invariant. Including 
the KQ75y7 extended temperature mask, we show the diagonal elements of the inverse noise 
covariance matrix in Fig. |5| 

Again, for the analysis, we generated 32 independent Monte Carlo chains with 140 000 
samples. After discarding the first 15 000 elements during burn-in, we applied the Gelman- 
Rubin convergence diagnostics to the chains and obtain a value of R = 1.5. The computed 
mean of (/nl) = 90 and the 1-a error of cr/j^^ = ^.re in agreement with the input values 
of the simulation. We show the constructed histogram on the right hand panel of Fig. [5} 
demonstrating the applicability of the algorithm to realistic data sets. 

6. Summary 

In this paper, we introduced an exact Bayesian approach to infer the level of non- 
Gaussianity of local type, /nl, from realistic CMB temperature maps. We derived conditional 
probabilities for the primordial perturbations given the data, P{^j,\d,9), and for /nl given 
the data and the perturbations, ^(/nlM, '^'l, We used Hamiltonian Monte Carlo sampling 
to draw valid realizations of $l from which we in turn sample /nl- After convergence these 
are samples from the full Bayesian posterior density of /nl given the data. 

For a direct comparison of the newly developed scheme to the conventional fast (bispec- 
trum) estimator, we used simulated Gaussian and non-Gaussian CMB maps superimposed 
by isotropic noise. Estimates of the error bars within the frequentist approach were derived 
from Monte Carlo simulations. As a result, we find consistent outcomes between the two 
approaches for the analyzed Gaussian map, in agreement with the fact that the fast estima- 
tor is optimal in the limit of vanishing non-Gaussianity. In the non-Gaussian case, however, 
the advantage of the exact Bayesian approach becomes important. Here, the uncertainty 
in /nl remains at the same level as for the Gaussian simulation, whereas the frequentist 
technique suffers from excess variance. Our results give the first example of an estimator 
(the "mean posterior estimator") that saturates the Cramer- Rao bound for /nl even if the 
signal is detectably non-Gaussian. 

Finally, we demonstrate the applicability of the newly developed method to a realistic 
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data set with spatially varying noise properties and partial sky coverage. Considering a 
WMAP-like noise covariance matrix and imposing the KQ75y7 extended temperature anal- 
ysis mask, we analyze a non-Gaussian simulation and recover the input value consistently. 

In the limit of undetectable non-Gaussianity, the Bayesian approach ought to yield 



the same information as the optimal bispectrum estimator (Babich 2005 Creminelli et al. 



2007). Even in that limit it is useful as a cross-check since it is implemented in a completely 
different way. Although being computationally expensive, we conclude that the method 
presented here is a viable tool to exactly infer the level of non-Gaussianity of local type from 
CMB radiation experiments within a Bayesian framework. 
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Fig. 1. — Properties of the maps analyzed. Left panel: Our non-Gaussian CMB signal 
simulation in dimensionless units. Right panel: The input signal (solid line) and noise 
(dashed line) power spectra. 




Fig. 2. — Analysis of the Gaussian simulation (/nl = 0). Left panel: We show the analysis 
of the Gaussian CMB map by means of the frequentist estimator. Plotted are the recovered 
value /nl = 4 {solid line) and the 2 — a error {dashed lines) as derived from Monte Carlo 
simulations, Right panel: The analysis of the same data set within a Bayesian 

framework constructs the full posterior distribution P{f]ssh\d, 9). We obtain a mean value of 
(/nl) = 3 and a standard deviation of af^^ — 15. 
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Fie;. 3. — Same as Fig. g but for the non-Gaussian simulation (/nl = 100). The results 
from a frequentist analysis are /nl = 97, crfj^'^ = 20. Using the Bayesian method, we obtain 
(/nl) = 99 and ctjj^^ = 15. For a significant detection of /nl, the bispectrum estimator 
shows excess variance, whereas the analysis on the basis of the exact Bayesian approach still 
provides tight error bounds. 
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Fig. 4. — Performance of the sampling algorithm. Left panel: We plot a random selection of 
three of the 32 /nl sampling chains that build up the histogram in Fig. |3| We discarded the 
first 10 000 samples during burn-in. Right panel: The autocorrelation function of a sampling 
chain. 
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Fig. 5. — Analysis in case of a realistic CMB experiment. Left panel: We show the diagonal 
elements of the inverse noise covariance matrix in dimensionless units adopted for the more 
realistic simulation. When expressed in real space basis, off-diagonal terms vanish exactly. 
Pixel within the KQ75y7 mask are set to zero, corresponding to assigning infinite variance to 
them. Right panel: The constructed posterior distribution P(/nlM, 0) of the simulated map. 
Obtaining (/nl) = 90 and af^^ = 17 for the mean and standard deviation, respectively, the 
input value (/nl = 100) gets consistently recovered. 



